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ABSTRACT 



This thesis examines acoustic transient discrimination and Time Difference Of 
Arrival (TDOA) estimation for the purposes of estimating the position of a submarine in 
a sonabuoy field. Transient discrimination, for this thesis, is the process of telling 
different transients apart. Two algorithms are evaluated. One method is based on higher 
order statistics while the other is based on signal subspace techniques. Extensive 
simulations using synthetic transients were conducted to establish the performance of 
each algorithm in terms of discrimination and TDOA estimation. It was found that the 
bispectral algorithm gave better TDOA estimation at low SNRs while the subspace 
algorithm gave better TDOA estimation at high SNRs. For discrimination, it was found 
that the subspace algorithm gave consistant false alarm rates at all SNRs while the false 
alarm rate for the bispectral algorithm grew with increasing SNR. 
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EXECUTIVE SUMMARY 



In this thesis we have developed and compared two algorithms, namely the 
bispectrum and subspace linear phase detectors. These algorithms were developed for 
the purposes of transient discrimination and Time-Difference-Of-Arrival (TDOA) 
estimation. They are to be used as part of a transient tool suite to aid in the estimation of 
a submarine’s position. Two performance measures were used to evaluate the 
algorithms, namely the probability of correct TDOA (Pt) and the probability of correct 
discrimination (Poi)- 

In general, it can be said that for TDOA, the bispectral linear phase detector gave 
better results at low SNRs while the subspace linear phase detector worked better at the 
higher SNRs. 

For discrimination, it was found that the bispectral discriminator gave higher Poi 
than the subspace discriminator. However, the probability of false discrimination (Pfa) of 
the bispectral discriminator increased at higher SNRs while the subspace discriminator 
gave a constant Pfa for all SNRs evaluated. For discrimination, the advantage of having a 
constant Pfa is desirable. Therefore the subspace discriminator is the best option even 
although it produced lower Pq; than the bispectral discriminator. It was also found that 
there are design trade-offs between processing speed and performance that need to be 
made. For the bispectral linear phase detector, this trade-off is in terms of threshold gain; 
for the subspace linear phase detector, this trade-off is in terms of correlation matrix size. 
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I. INTRODUCTION 



The work of this thesis forms part of the ongoing effort to automate the detection, 
discrimination and Time-Difference-Of- Arrival (TDOA) estimation of transient signals. 
At present, transient signals are located and processed manually, making it a very time 
consuming procedure. It is therefore more desirable to have an automated system that 
can deliver the same or better results than the present manual approach. 

For this thesis, it is assumed that the detection of the transient has already taken 
place; therefore, the objective of the thesis is to look at discrimination and TDOA 
estimation only. Two algorithms, namely the bispectral linear phase detector and the 
subspace linear phase detector, are developed for this purpose. 

A. THE OPERATIONAL MISSION 

A typical mission considered for this thesis starts off with an anti-submarine 
warfare (ASW) aircraft laying a sonabuoy field in the vicinity or on the predicted path of 
a submarine. Usually each sonabuoy in the field consists of calibrated omni-directional 
passive acoustic sensors. The ASW aircraft collects acoustic data on a target submarine 
as it passes through the sonabuoy field. 

Figure 1 shows a typical V-shaped sonabuoy field that can be used in this type of 
mission. When this shape is used, the target must pass through the apex of the V of the 
field (as close as possible) in order to obtain the most accurate results. The data received 
at the sonabuoys is transmitted to the aircraft, which in turn records the data on tape for 
mission post processing. 
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The data collected on the tapes is processed at an onshore location in order to estimate the 
acoustic Sound Pressure Levels (SPL) of tonals, broadband signals and transients that are 
emitted from the submarine. 

To obtain precise SPL of a submarine, an accurate track is required so that the 
pressure levels received at the sonabuoys can be projected back to estimate the pressure 
levels at the target. For best results, the position of each sonabuoy must be known, the 
target must be on a course that passes through the middle of the field, and an accurate 
estimation of TDOA of acoustic signals must be made. 

B. DATA PROCESSING 

The data collected from the sonabuoys and stored on the tape is processed by first 
generating a baseband sonogram. A sonogram is a frequency versus time plot of the 
acoustic pressure levels at a sonabuoy. From these sonograms an analyst can identify and 
highlight contact events of interest, which may include signals, such as narrowband 
tonals and transients. The events identified are then processed to obtain an estimate of 
the target’s track. Once the best estimate of the target track has been made, the analyst 
can conduct SPL analysis on all types of previously identified contact signals. The final 
desired results are the SPL signature characteristics for a given target with respect to 
aspect angle. 
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At present, transient analysis is a mcinually intensive operation. Typically, the 
analyst marks transient events of interest from the sonogram display of a sonabuoy and 
then searches the sonograms of the other sonabuoys for the same transient. This can be a 
very time consuming process since in a single mission there can be hundreds of transients 
within a single sonogram. These must be matched up to those on the sonograms from the 
other sonabuoys in the sonabuoy field. 

C. THESIS GOAL 

The goal of this thesis is to develop and evaluate two algorithms that can be used 
for transient discrimination and Time-Difference-Of- Arrival (TDOA) estimation. The 
two algorithms are the bispectral linear phase detector and the subspace linear phase 
detector. 

D. THESIS OUTLINE 

The remainder of this thesis is organized as follows. Chapter II discusses 
transients and transient processing and introduces the problems of interest, namely 
TDOA estimation and transient discrimination. Chapter IQ presents the relevant signal 
models and the analysis of the measured signals in terms of second- and third-order 
moments. This chapter details the formulation of the bispectral linear phase detector and 
discriminator. Chapter IV discusses signal subspace techniques and their application to 
the problem of TDOA estimation and transient discrimination. The outcome of this 
chapter is the formulation of the subspace TDOA estimator and discriminator. Chapter V 
examines the results achieved by both the bispectral and subspace TDOA estimators and 
discriminators. Extensive simulations are conducted using synthetic transients to produce 
performance curves and to study the utility of each technique in discrimination and 
TDOA of transients. Chapter VI provides conclusions and recommendations for future 
work based on the results presented. Appendix A and Appendix B show the 
discrimination results for the different combinations of signals not shown in Chapter V. 
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II. TRANSIENT PROCESSING 



A. TRANSIENT SIGNALS 

Transients of interest here are short duration wideband acoustic signals. 
Accordingly transients can be of varying shapes and lengths, lasting anywhere from a few 
microseconds to a couple of seconds. Typical examples of these signals are a wrench 
falling on a metal deck or the sounds of a biological life form, such as a whale. Due to 
the diversity of these signals, there is limited a priori information that can be used to aid 
in their detection. When transients are imbedded in noise, their detection becomes very 
difficult because the energy in the noise frequently dominates the energy in the transient 
over the interval of interest. Nevertheless, in conjunction with other measurements or by 
themselves, transients can be used to estimate the track of the target submarine by 
estimating the TDOA of a transient between two sonabuoys. 

Tracking of a submarine using transients requires transient detection, 
discrimination and TDOA estimation. Prior to discrimination, a transient must be 
detected in a noisy background on the sonagram of a single sonabuoy. This thesis does 
not address the transient detection problem and therefore assumes that a transient has 
been detected on the sonogram of a sonabuoy. 

After the transient has been detected, the same transient must be found on the 
sonograms of the other sonabuoys in the sonabuoy field. This is difficult since the 
detected transient may not exist on all sonabuoys within the field. That is, the sonograms 
of the other sonabuoys could contain only noise or even a different transient in the 
window of interest. The first case (i.e., where the sonogram of the other sonabuoys 
contain only noise) is a common scenario. This happens because transients can be very 
localized and thus do not appear on the sonograms of each sonabuoy. The second case 
(i.e., where a different transient exists on the sonogram of the second sonabuoy) is also 
common due to the complex environment of the mission, where multiple transients from 
different sources can exist on the data records of the other sonabuoys in the field. It is 
therefore important to discriminate between different transients and to know if there is no 
transient present. If the transients arriving at two or more different sonabuoys are the 



5 



same, then TDOA estimates can be used for target localization. In the following we first 
discuss TDOA estimation and then transient discrimination. 



B. TIME DIFFERENCE OF ARRIVAL (TDOA) ESTIMATION 

TDOA estimation is by no means simple and therefore many authors have written 
about it [Ref 1]. It is, however, the primary means of determining range to the target in 
passive detection, since TDOA information from multiple sonabuoys and the geometry of 
the sonabuoy field can be used to determine the position of the submarine at any 
particular instance in time. For any given value of TDOA between two sonabuoys, the 
locus of possible target positions is a hyperbolic curve. If more than one curve can be 
drawn, i.e., if TDOA values between three or more sonabuoys can be calculated, the 
intersection of the curves determines the position of the target. 

Consider a simple two-buoy model as shown in Figure 2. The sonabuoys are 
symmetrically positioned about the origin on the x-axis as shown in Figure 2, with the 
target located at T. 



Figure 2. Two Symmetrically Positioned Sonabuoys. 

For this geometry the time delay t between sonabuoy 1 and sonabuoy 2 is given 




Buoy 2 (-xo,0) Buoy 1 (xo,0) 



by 



I ‘■2 

C C 




( 2 . 1 ) 



{r-cf =(V(^o +y^ +y^y 
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where c is the speed of sound, D,^ and D,^ are the distances between the target and 

sonabuoy 1 and sonabuoy 2, respectively. After some algebraic manipulation’ of Eq. 2.1, 
the following expression is obtained [Ref 1]: 



£1-Z^ = i 
a' 6' 



( 2 . 2 ) 



where 



T.C 

a = — 
2 




Equation 2.2 describes a hyperbola, and a set of hyperbolas can be drawn for different 
TDOA’s between sonabuoy 1 and 2 as shown in Figure 3. 




Figure 3. Hyperbolic curves for TDOA between two symmetrically positioned 
sonabuoys. Sonabuoy positions are indicated by the dots. 

C. TRANSIENT DISCRIMINATION 

The objective of discrimination between transients is to be able to tell different 
transients apart. This is different from transient classification or transient 
characterization, which attempts to identify the source of the transient and seeks to group 
similar transients together. 



' To obtain a hyperbolic expression, Eq 2.1 must be squared, and all the cross terms must be left on the 
right hand side and all other terms on the left hand side. Both sides must now be squared again with like 
terms collected. 
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Discrimination is difficult since it relies on detecting differences between the two 
received signals. These differences can be either in magnitude or phase or both. The 
problem is further complicated by the presence of noise and the fact that (due to 
differences in propagation path characteristics) the same transient arriving at two 
different locations may not look and sound the same. Achieving good discrimination at 
low SNR values is a challenging task. 

Figure 4 shows some typical cases in which transient discrimination has to take 
place. The cases shown in this figure are a localized transient such as an expanding 
bubble on a sonabuoy, a directional transient that is emitted from the hull of the 
submarine and is only received at some of the sonabuoys in the sonabuoy field, and an 
omnidirectional transient that is received at all sonabuoys in the field. 




Local Transient 



o 



o 




o 




Figure 4. Discrimination Cases. 



Considering these cases, the TDOA cannot be estimated for all combinations of 
two sonabuoys in the field. Further, if a transient arriving at one sonabuoy is mistaken to 
be the same as the transient arriving at another sonabuoy, the resulting TDOA estimate 
will produce an erroneous target position. This in turn may lead to errors in tracking and 
ultimately possible erroneous SPL calculations. 
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III. MOMENT-BASED SIGNAL PROCESSING 



Most of the traditional work on transient processing and TDOA estimation is 
based on the use of second order statistics and classical (i.e., Fourier-based) methods of 
spectrum estimation [Ref 2], [Ref 3], [Ref 4], [Ref 5], [Ref 6]. Recently, new work has 
appeared using techniques involving higher order moments of signals and higher order 
spectra, such as the bispectrum [Ref 7], 

The motivation for using the bispectrum in transient discrimination and TDOA 
estimation is twofold. First, the higher order spectra suppress Gaussian noise processes 
of unknown spectral characteristics. Secondly, these spectra, unlike the usual power 
density spectrum, preserve phase information [Ref 7]. In the last few years there has 
been a considerable amount of new research done in using the bispectrum and 
trispectrum for transient detection, time delay estimation and classification of signals 
[Ref7;p 313], [Ref 8], [Ref 9], [Ref 10], [Ref 11], [Ref 12]. 

This chapter defines the signal model used in this thesis and the analysis of these 
signals based on second and third moments. Finally, a set of algorithms for 
discrimination and TDOA estimation using third order moments is discussed. 

A. SIGNAL MODEL 

Before any further analysis can be done it is necessary to develop a signal model, 
which can be used for the analysis of transients arriving at two sonabuoys. As mentioned 
previously, there are three cases to be considered. In the first case, the transient arriving 
at the second sonabuoy is the same as the transient arriving at the first sonabuoy. In the 
second case, the two arriving transient signals are different while in the third case only 
noise is present at the second sonabuoy. It is also assumed that the sonabuoys used are 
omni-directional and that their separation in the sonabuoy field is large enough so that 
noise at the first sonabuoy is uncorrelated with noise at the second sonabuoy in both 
space and time. 

The following simple model can be used to represent a single transient arriving at 
two different sonabuoys: 



9 



(3. 2) 



xXk) = s{k) + T}Xk) 

Xj (A:) = As{k - 1) + 7 , (A:), 

where ^, (A:) is the noise-embedded signal arriving at ^onabuoy i, ^(A:) is the transient 
signal itself, and Tj.{k) is white gaussian noise. Note that the signal at sonabuoy 2 is 

subject to a relative attenuation A and delay L with respect to the signal at sonabuoy 1. 
The frequency domain expression for these signals is 

xXco) = S{(o)+ N,{co) 

{a>) = + N2 (<y), 

where the uppercase letters represent Fourier transforms of the respective signals in Eq. 
3.1. 

If the transient arriving at sonabuoy 2 is different from the transient arriving at 
sonabuoy 1, the received signals are given by: 

x,{k) = s{k)+Tj^{k) 

X2{k) = r{k-L)+T]2{k), 

where r(A:) is the other transient arriving at sonabuoy 2 with time delay L relative to s{k). 
The corresponding expressions in the frequency domain are 

^,(<y) = 5(6))+iV,(<y) 

X2{co) = R{(o)e-^‘^ +N2 {co). 

If there is no transient present at the second sonabuoy, the two received signals 
are 

xXk) = s{k)+n,{k) 

X2[k) = rh{k)- 

In this case the signal received at sonabuoy 2 consists entirely of noise and the frequency 
domain representation of the two signals is thus: 

X^ (co) = ^(^y) -I- Nj (co) 

X,(<y) = A^,((y). 



(3. 3) 



(3.4) 



(3. 5) 



(3. 6) 



B. MOMENTS AND CUMULANTS OF A RANDOM PROCESS 

For the purposes of this thesis, the following notation and definitions are used. 
The n"‘ moment of a real stationary random process is [Ref 7:p. 15] 
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E{x{k)x{k + r,)..jc{k + r„_^)\ (3. 7) 

where E denotes statistical expectation and r,. represents the lag. Note that w^(r) is 
the ordinary autocorrelation function . 

In the analysis of signals using higher order statistics, cumulants rather than 
moments are generally used. Cumulants of order n are defined by certain linear 
combinations of products of moments of order n and lower [Ref 7:p. 15]. Cumulants of 
Gaussian random processes have the distinction that they are all zero for order n greater 
than 2. Thus signals imbedded in additive gaussian noise theoretically appear naked 
when subjected to analysis using higher order cumulants. The order cumulant is 
denoted by 

,v,,)=CuOT[x(k)x(^ + r,)..Jc(A: + r„.,)] (3. 8) 

For a more detailed definition of cumulants the reader is referred to [Ref 7:p. 9]. 

The following relationships exist between moments and cumulants for stationary 
random processes [Ref 7:p. 9]: 

a. First Order. 



c;=n,;=£{4*)) (3.9) 

b. Second Order. 

c'(r,) = w;(r,)-(m])' (3.10) 



c. Third Order 

c\ (r, ,T^) = m] (r, ,T,)-(m\) [m] (r, ) + m] {r, ) + (r^ - r, )] + l(m\ J 

(3.11) 

Using these relationships, zero mean, white Gaussian noise has the characteristics listed 
in Table 1. 
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n 




c" 


1 


0 


0 


2 






3 


0 


0 


4 


3a^ 


0 



Table 1. Moments and cumulants for white gaussian noise for r, = Tj = Tj = 0. 



C. N-TH ORDER “MOMENTS” OF A DETERMINISTIC SIGNAL 



In this thesis, we consider transients to be deterministic signals. Consequently, 
concepts such as statistical moments are not defined. However certain operations in the 
time domain, which are analogous to estimating moments for realizations of a random 
process, are still useful for deterministic signals. Some authors (e.g., Nikias and 
Petropulu [Ref 7]) have referred to these operations as computing “moments and higher 
order spectra for deterministic signals.” Since some of the techniques we have adapted 
are due to authors using this concept, we will adopt this concept here as well. 

In general the n'* order moment for an energy signal, a:(^), is defined as [Ref 7:p. 
78] 



00 

»i"(r,,....r„.,)= S x{k)x{k + T,}..x{k + T„_^), 

k=-co 



and for a power signal [Ref 7:p. 100] as 

J J + N-l 

I xik)x{k + T^)....x{k + T„_^), 

1\ k=J 



(3. 12) 



(3. 13) 



where J is an arbitrary starting point of the summation and N is the period of the signal. 
From these expressions, moments can be considered as a measure of the degree of 
similarity between a signal and delayed or advanced replicas of itself. The n'* order cross 
moment for n energy signals is defined as 

(^-I v-,r„.,) = I xXk)x^{k + r,)...x„(^ + r„.,) (3. 14) 
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The n‘^ order moment spectrum is the Fourier transform of Eq. 3.12 and can be 
expressed as Fourier transforms of the signal. For energy signals, this is given by [Ref 
7:p.85] 

= X{co,)..X{co„_^)X\q}, + (3. 15) 

which can be written in terms of magnitude and phase as 

M; (®, )| H X(co , ) I .. I X(co „., ) II X(0), +.. + co „., ) | 

(3. 16) 

4^; (<y, ) = </>M+... + ), 

where ,...<y„_,)| is the magnitude term and 'F^"(iy,,...o)„_,) is the phase term. 

Orders n=2, 3 and 4 are important special cases of moments. In the frequency domain 
these orders have been termed Power Spectral Density (PSD), Bispectmm, and 
Trispectrum, respectively. 

D. SECOND-ORDER MOMENTS 

1. Definitions of Moments and Spectra 

Using Eq 3.12, the autocorrelation of a deterministic signal, x(A:), is written as 

E x{k)x{k + t\ (3. 17) 

k=-oo 

From Eq 3.16, the magnitude and phase components of the PSD are given by 

•v;{a)=0. 

From these expressions, it can be seen that the PSD is an even function with no phase 
information. Similarly, the cross-correlation of two deterministic signals is 

^ xXk)x^{k + T). (3. 19) 

*=-oO 

The corresponding cross-spectrum is the Fourier transform of Eq 3.19: 

(3.20) 

which can be written in terms of magnitude and phase 



(3. 18) 
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(3.21) 



2. Second Order Moments of the Received Signals 

Since cross-correlation is used extensively throughout the thesis, it is important to 
determine the cross-correlations for the three signal cases presented earlier. Before doing 
this, let us first investigate the cross-correlation between the two noise sources 7 , (A:) and 
at the two sonabuoys of interest. The cross-correlation and the corresponding 
cross-spectrum are given by: 

k=-oo (J. ZZ) 

K„ (®) = sK,. (r))= N, («K (ffl) 

The expectation of this term is zero since the two white noise sources are uncorrelated. 
However, the term, , itself is not zero. Applying this to the case where the same 
transient arrives at the two sonabuoys, the cross spectrum is 

= . (3.23) 

where are the cross terms. From this equation, it can be seen that the time delay L 

is imbedded in the linear phase term . This linear phase term is therefore important 
in finding the TDOA between the two signals. 

For the case where different transients arrive at the two sonabuoys, the cross- 
spectrum is 

= (3. 24) 

For this case it can be seen that there is the same linear phase term as that of Eq 

3.24. Therefore the same TDOA result would be observed for Eq 3.24 and 3.25 even 
though the transients are different. The linear phase can therefore not be used as a means 
to discriminate between two transients. In the last case, where there is only noise present, 
the cross-correlation is 

(<») = <„ (4 ( 3 . 25 ) 
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For this case, it is expected that the phase will be random, and therefore no linear phase 
term will exist. 

E. THIRD-ORDER MOMENTS - BICORRELATION AND BISPECTRUM 
1. Definition of Moments and Spectra 

The third-order moment of a signal is called the bicorrelation. From Eq 3.12, the 
auto-bicorrelation of a deterministic real signal is defined as 

(^I > ^2 ) = 2 + ^1 + ^2 ) (3.26) 

k--oo 

and fi'om Eq 3.14 the cross bicorrelation is 

"h', X,., (^1 , ^2 ) = 5 )-^3 (^ + ^2 ) (3- 27) 

For this thesis, where there are only two data streams, the cross-bicorrelation will consist 

of only two signals, Xj (A:) and (k) . The cross-bicorrelation functions of interest are 

00 

2: X,(k)x,(/c + rJx,(k + T,) 

k^-oo 

and the terms > '^xjx.x, > which are defined in an analogous way. 

The auto-bispectrum is 

Ml {co^ ,0)2)= X{co^ )x{co 2 )X'{q}j + CO2 ), 

or in terms of magnitude and phase 

|m> I = p(a), + <», )| 

'f'i (<», . ) = «>, (®, )+ «>. (‘“2 ) - «>. ('f 1 + “ai ) 

As can be seen, by comparing Eq 3.18 and Eq 3.30, a distinct difference between the 
second order moment spectrum and the bispectrum is that the bispectrum contains phase 
information. The cross-bispectrum is 

+®2>. p. 31) 

or in terms of magnitude and phase 



(3. 28) 



(3. 29) 



(3. 30) 
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W„io,_'\ = \X,_{a,]\X, (o, JA", (o, + a>, )| 

w, (®l . ) = !^., (‘»l ) + «>,. (®! )-«>,, (<i>l + <»i ) 



(3. 32) 



2. Third Order Moments of the Received Signals 

Both second order moments and third order moments are used extensively 
throughout this thesis. It is therefore important to develop the third order moments for 
the three signal cases. To provide a motivation for the algorithms to be used, let us 
consider a situation in which the noise is identically zero. In reality the noise terms 
N^(co) and A^ 2 (^) zero, although the expected value of their higher order 

moments defined earlier vanish when the noise is Gaussian. 

For the case where the same transient arrives at both sonabuoys, the signals in the 
frequency domain are given by Eq. 3.2. Using Eq. 3.32, the bispectrum of these signals 
is given by 

("l ’ ^^2 I = F(^^2 +(^ 2 } 

'i' (". ,co^)=^sM-o)L + (f>s {0)2 ) - <f>s ("1 + " 2 ). 
where (j>^ is the phase of the signal 5(<y). Equation 3.33 is the expression for the cross- 
bispectrum between the data streams at sonabuoy 1 and sonabuoy 2 in the absence of 
noise. The signal phase terms, ^j(^i)+^s(^ 2 )“^s(*^i ■'■‘^ 2 )’ ^e eliminated by 
making use of the auto-bispectrum at sonabuoy 1, which is given by 



(3. 34) 



1^ i (^^1 ’^2 )| = +^2} 

(^^1 ^^2) = <PsM+ <Ps ("2 )- <!>S (^1 +<^ 2 )- 

This expression is then used to normalize the cross-bispectrum and results in the identity: 

,(^^ 1 .^ 2 ) 












(3. 35) 



From Eq 3.35 it can be seen that by normalizing the bispectrum all the phase terms are 
cancelled except for the linear phase term, . 
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For the case of different transient arrivals, the cross-bispectrum can be obtained 
from Eq. 3.4 and Eq. 3.32 as follows 



)| = +(^ 2 } 

'i'.M {o>x ,co^) = (I>rM-o3L + </>s {co, ) - {(O, + ), 

where ^^{co) is the phase of the signal R{co). In this case the normalized cross- 
bispectrum takes the form 



T(fn rn W ^ _ l^(^i )| 

i\a>^,W2)- 3 / \ “ lo/' \ ^ ■ (4. J/) 

M^\co„(02) |*S(tyij| 

Note that, unlike the first case (see Eq. 3.35), the linear phase term caimot be separated 
from the other phase terms. In both cases however, the two-dimensional bispectrum is 
reduced to a function of a single variable d>, . 



F. SIGNAL PROCESSING ALGORITHMS 



1. Bispectrum Linear Phase Detector 

The analysis of the previous section shows that when the same signal is arriving 
at both sonabuoys, a linear phase term is present in the normalized function l(o}j,£ 02 )- 
The time delay L can now be extracted using the following ad hoc method developed for 
time delay estimation [Ref 16]. 

To estimate the delay L, a third-order “hologram” transformation is required. 
This is defined to be [Ref 7:p. 324] 

n n 

T{t) = I \l{co„co2 y^‘^'^do),dco2. (3. 38) 

-n~K 

Since I reduces the bispectrum to one dimension, the third-order hologram is a one- 
dimensional Fourier transform over the dimension containing the linear phase term, 
followed by an integration over the second dimension. Note that when the same signal is 
present at both sonabuoys, Eq. 3.38 takes the form 
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(3. 39) 



T{^)= J ^^dco^dco^. 

-n-K 

Since all the elements of this second dimension are in phase they add up constructively, 
giving a strong peak at r = Z,. Since we are working in the discrete-time domain the 
third-order hologram can be rewritten as 

(3. 40) 

<y,=0<W2=0 



The absolute value of J(r) will display a strong peak at the location of the of the time 
delay between the two signals of sonabuoy 1 and sonabuoy 2. 

For the case where the two transient arrivals are different, the third-order 
hologram takes the form: 









y}<^\ 



{L- 






(^t )) 



dco^dco^- 



(3.41) 



In this case, the third order hologram contains extra phase terms <f>R and (p$, which will 



either add in phase or out of phase. Thus, in general, T(r) will not exhibit a strong peak. 



2. Bispectrum Discriminator 

The bispectrum linear phase detector can also be used as a discriminator by 
applying a simple threshold technique to the third order hologram. The threshold 
technique is applied by first noticing, from Eq 3.41, that if the transients are different the 



hologram contains phase terms ^-coL and magnitudes 






On the other hand, if 



the transient signals are the same, the hologram contains only the linear phase term coL 
(Eq 3.39). The phase and magnitude terms of Eq 3.41 will add constructively or 
destructively to produce peaks and valleys to the bispectral hologram. If the signals are 
the same, there will be only one peak in the hologram, and that peak will be at the delay 
L. 

In summary, one can expect that if the SNR is sufficiently large and the signals 
are the same, a peak at delay L will dominate the hologram. On the other hand if the 
signals are different, there will be a maximum value at delay L as well as other extrema at 
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delays other than L. These other extrema will be larger than the noise and can therefore 
be used to discriminate. A discrimination algorithm is therefore proposed as follows: 

Step 1: 

Estimate T{t) 

Step 2: 

Find the maximum value of the hologram 



2 = argmax T{t) 



(3. 42) 



(Q is the magnitude of the hologram at the estimated TDOA). 

Step 3: 

If Tq is the value of r that produces the maximum in step 2 then compute 




T 



(3. 43) 



o 



Step 4: 

Define a threshold using Eqs. 3.42 and 3.43 as 




(3. 44) 



Step 5: 

Compare Qj. to the difference Q —Q 2 - If 



Qt>Q-Q: 

the signals are considered to be different; if 

Qt<Q-Q: 

the two signals are considered to be the same. 



(3. 45) 



(3. 46) 



19 



THIS PAGE INTENTIONALLY LEFT BLANK 



20 



IV. SUBSPACE-BASED SIGNAL PROCESSING 



Signal subspace methods have been used in a variety of problems for estimating 
spectral lines and direction of arrival in sensor arrays. A good introduction to these 
methods for spectrum estimation can be found in [Ref 17]. Recently, the use of subspace 
methods has been explored for estimating TDOA [Ref 20], [Ref 21]. This thesis uses the 
same approach as [Ref 20] for TDOA estimation and further adapts the subspace method 
for transient discrimination. 

The discussion begins by defining the subspace model for the TDOA problem and 
then moves on to the MUSIC method, which was used for delay estimation and 
discrimination. The chapter finishes with a discussion of how the MUSIC method was 
adapted and used for discrimination 

A. SUBSPACE 

To understand the principle of subspace methods and their application to the topic 
of this thesis, consider the case where the same signal is arriving at two separate 
sonabuoys (see Eq. 3.1). The frequency domain representation of Eq. 3.1 is given by Eq 
3.2. From Eq. 3.2, it can be seen that the linear phase term contains the delay L 
between the transient arrivals at sonabuoy 1 and sonabuoy 2. The cross-spectrum for the 
two received signals, given by Eq 3.23, is repeated here for convenience 

{<•>)■ ( 4 . 1 ) 

It can be argued [see Ref 20] that when Eq 4.1 is sampled at equally- spaced values in 
frequency, the resulting data sequence, yik) = satisfies the conditions 

necessary to apply a signal subspace model [Ref 17]. Specifically, an AxA correlation 
matrix is estimated for the datay(^)> the corresponding A-dimensional vector space is 
divided into signal and noise subspaces. The approach followed to estimate TDOA using 
subspace methods is to project a steering vector of the form 
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1 



e 






d(/)= . 



(4. 2) 



e 



JUN-l) 



onto the noise subspace and plot the result as the linear phase, /, is varied. When the 
parameter / is equal to the true delay L between the sonabuoys, the projection onto the 
noise space is zero. In subspace algorithms, such as MUSIC, this null projection can be 
used to estimate TDOA. 

Unfortunately, the cross-spectrum for the case where different transients are 
arriving at sonabuoys 1 and 2 also contains a linear phase term (see Eq. 3.24). Therefore, 
the basic subspace techniques are not useful for transient discrimination without suitable 
modification. 

B. MUSIC 

In the application of subspace techniques to linear phase detection between two 
transients, the MUSIC (Multiple Signal Classification) method developed by Schmidt 
[Ref 18] is used. A brief outline of the method is provided here. 

The MUSIC algorithm is implemented in the frequency domain by first obtaining 
the cross-spectrum of the two signals received at the sonabuoys. Since is 

formed using a DFT, this data is available at samples , k= 0,1,2 ...Ndft-L where Ndft 
is the size of the DFT. From this frequency domain data, a data matrix^ is formed, 
and the correlation matrix is estimated as 




(4. 3) 



This can be decomposed into an orthonormal eigenvector matrix 



and a diagonal eigenvalue matrix 



(4. 4) 



^ See [Ref 17] for various methods of forming a data matrix 
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A = 



0 



(4. 5) 



The columns of Ej/o and E„o„e form an orthogonal basis set and define the signal and 
noise subspaces, respectively. It is therefore possible to write as 

R, = EAE*^ = E_A„. eX + E„,;,,A„_E'L- (4. 6) 

X sig sig sig noise noise noise v / 



The columns of the eigenvector matrix can be used to form projection matrices for the 
signal and noise subspaces of the form [Ref 17:p. 623]: 



=E„,, E’L. 



(4. 7) 



When the projection matrix P is multiplied by a vector d, the result d = Pd is the 
projection of d into the corresponding subspace (see Figure 5). 




Figure 5. Projection of vector d onto noise subspace by MUSIC. 



As discussed previously, subspace methods make use of the fact that the 
projection of the steering vector d(/) , given by Eq. 4.2 onto the noise subspace is zero. 
The MUSIC algorithm [Ref 17:p. 627], in particular, evaluates the quantity 



P MU 






d-(/Ee,erd(/) 

i=2 



(4. 8) 
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where is the MUSIC indicator function for d(/)and e 2 -..e^ are the eigenvetors 
spanning the noise subspace (the eigenvector e, spans the signal-subspace which is one- 

A 

dimensional) The value of / where exhibits a sharp peak determines the TDOA 
between the signals. 

C. SUBSPACE DISCRIMINATOR 

A subspace discriminator was developed by comparing projections onto multiple 
subspaces using the MUSIC algorithm. The following considerations are proposed to 
motivate this discriminator. 

Since the subspace method in this thesis is based on the cross-spectrum between 
two sonabuoys (see Eq. 4.1) it is important to investigate the differences between the 
cross-spectrum for the case where the same transient arrives at two separate buoys and 
that for the case where different transients arrive at two separate buoys. This will be 
helpful in understanding the subspace discriminator. From Eq. 3.24, it can be seen that if 
the transients are different, the cross-spectrum contains phase terms ^g{p))-(j)^{co)-coL 

and magnitudes |p(<y,)||5(<y,)| . However, if the transient signals are the same, the cross- 
spectrum contains only the linear phase term (oL (Eq. 3.23). Now, if the subspaces 
formed using Eq. 3.23 are compared to those formed using Eq. 3.24 the two subspaces 
would typically be rotated with respect to each other as illustrated in Figure 6. The 
implication of the rotation of the subspaces due to the phase terms of Eq. 3.24 is that 
when the parameter / of the steering vector, d(/), is equal to the delay L between the 
sonabuoys, the projection onto the noise subspace will not be zero as in the case when the 
same transient arrives at the two sonabuoys. By comparing projections, it would be 
found that the differences in magnitude of the projections would be large if the transients 
were the same. On the other hand if the transients were different, the differences in 
magnitude between the projections would be small. These differences in magnitude are 
used to discrirrmate between transients. 

This method is developed further by noticing that the noise projection matrix, 
^noise, can be written as: 
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(4-9) 

1=2 

where I is an NxN identity matrix.Let us further define P* as the projection operator for 
the subspace formed by leaving out the k''’ basis vector. Thus 

P. =Se,e:"=I-e.e7. (4.10) 

1=1 

i*k 



Further we define the following set of additional indicator functions 

1 1 






d-"(/)P,d(;) d-"(/Xl-e.e7]d{/)’ 

The proposed sub space-based discrimination algorithm is as follows: 



k = 2,...,N. (4.11) 



Step 1 

Estimate the correlation matrix Rx and find the eigenvectors e,. 
Step 2 

Compute the maximum estimate 



Q = argmax 

k.l 



f 

P 

J 



Step 3 

If ko is the value of k that produces the maximum in step 2 then compute 

Q 2 = argmax 

k.l 

k^kQ 

Step 4 

Compare the difference Q -Q^- If 

e-^2<l 

the signals are considered to be different. If 

Q-Q 2 <1 

the signals are considered to be the same. 




(4. 12) 



(4. 13) 



(4. 14) 
(4. 15) 
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Figure 6. Relative rotation of subspace due to phase terms of Eq. 3.24. 
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V. SIMULATION RESULTS 



A. SIMULATION CONDITIONS 

1. Synthetic Transients 

In this thesis a set of four synthetic transients was used to evaluate each of the 
algorithms described in the previous chapters. The four transients were a CW sinusoidal 
pulse, an exponentially decaying sinusoidal pulse, a linear phase modulated (LFM) pulse 
and a simulated finback whale transient. 

Each transient can be described as an amplitude- and phase-modulated sinusoid, 
given by the general expression [Ref 19:p. 202] 

p{t) = a{t)cos6{t), (5.1) 

where a(t) is the time-varying amplitude or “envelope” of the signal and 6(t) is the time- 
varying angle. The angle 6(t) is of the form 6{t) - + y{i) so that 

p{t) = a{t) cos{o)j + y{t)), (5. 2) 

where fc is the center frequency and y(t) is the phase modulation. Each of these transients 
is described in more detail below. 

For the CW pulse, the envelope and phase modulation are given by: 

fl, 0<t<T 

a{t) = \ (5. 3) 

[0, otherwise 

y(0 = 0 (5.4) 

where T is the pulse length. A plot of this transient is given in Figure 7 for T = AO ms 
and /, =50 Hz. 
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Figure 7. CW Pulse. 



The exponentially decaying sinusoidal transient is described by the following 
envelope and phase modulation terms 

fe’", Q<t<T 



a{t) = ■ 



0 , 



otherwise 



(5. 5) 



r(0 = 0, (5.6) 

where T is the pulse duration and a is the attenuation constant. This transient is shown in 



Figure 8 for a = 200 and T = 40ms . 
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Figure 8. Exponentially Decaying Sinusoidal Transient. 



For the Linear Phase Modulated (LFM) transient, a Gaussian envelope of the 



form 
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(5.7) 



a{t) = 



i-iV 

e , 0 < ; 



0 , 



<T 
otherwise 



was used. The phase is a quadratic function of time: 

Ik • Isf 



yit) = ■ 



IT 



(5. 8) 



where A/" is the desired change in instantaneous frequency over the interval T. A plot of 

T 

this transient for T = AOms ,< 7 ^ = — and Af = 1007/2 is shown in Figure 9. 
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Figure 9. LFM Pulse. 



The finback whale transient is the most complicated of all the transients that were 
synthesized. This transient is modeled as follows [Ref 8]: 

\3IT)-t, 0<t<7’/3 

1, r/3<r<2r/3 

a{t) = . (5.9) 

3-(3/7’)-f, 2T!3<t<T 

0, otherwise 

y(0 = 2;r-23.e'”('*‘''''‘‘^' (5.10) 



with fc = 0. This transient is plotted in Figure 10 for T = 1 5 . The instantaneous 
frequency decreases nonlinearly from 23 Hz to 18 Hz with the most rapid decrease 
occurring initially and the rate of decrease becoming smaller with time. 
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Figure 10. Finback Whale Transient. 

2. Signal-to-Noise Ratio 

During the evaluation of the algorithms, the synthetic signals described above 
were subject to additive white Gaussian noise. The signal-to-noise ratio (SNR) is defined 
as follows. Let represent the average power in the signal: 

Kcnsien, = ^ £k(«) T ' (5- 1 1) 

n=0 

Note that the transient signal power is normalized by the length of the transient Ns and 
not the entire observation time. This is done to prevent the SNR from changing with 
changing observation time. The SNR in dB is then defined as 

SNR = lOlog.o = lOlog.o dB (5.12) 

P CT^ 

Noise noise 

where is the noise variance used in the generation of the Gaussian white noise. 

B. RESULTS 

To evaluate the algorithms, two sets of experiments were used, one to perform 
TDOA estimation and the other to perform discrimination. For the first set of 
experiments, the same transient arrives at both sonabuoys. For this case all the transients 
were evaluated in turn, over a range of SNR values. For the second set of experiments. 
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all combinations of transient arrivals at sonabuoys 1 and sonabuoys 2 were tested for 
discrimination. 

In all of the experiments, the transients were set to the following lengths: 20 ms 
for the CW sinusoidal pulse, 40 ms for the exponential decaying sinusoid and LFM pulse, 
and 1 s for the whale transient. The delay, L, between the transients was always kept at 
40 samples, which is equivalent to 200ms. A sampling rate of four times the required 
Nyquist frequency was used for each transient. This satisfies the requirement for the 
bispectrum, which requires a sampling rate of three times the maximum frequency [Ref 
8]. The total length of each signal used was 256 samples or 1.28 s. The parameters are 
summarized in Table 2. 





/.(Hz) 


T 


L (delay) 


A/ 


a 


CW Pulse 


50 


20ms 


200ms 






Exponential 
Decaying Sinusoid 


50 


40ms 


200ms 




200 


LFM 


50 


40ms 


200ms 


IkHz 




Whale 




Is 


200ms 







Table 2. Transient parameters used in experiments. 



1. TDO A Estimation 

The TDOA results will be discussed as follows. First, an example of the basic 
results using the exponentially decaying sinusoid transient will be given for both the 
subspace method and the bispectrum linear phase detector. After this, certain 
“probability” curves for TDOA will be defined and presented for each transient starting 
with the subspace method and ending with the bispectrum results. 

The subspace method uses the MUSIC algorithm with a correlation matrix (Eq. 

4.3) of size N = 60. Figure 1 1 shows a plot of the function Pmu (see Eq. 4.8) as a 
function of / using a SNR of 12 dB. Figure 12 is a plot of the results for a SNR of 5 dB. 
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Figure 11. Subspace TDOA estimation for an exponentially decaying sinusoid 
using an SNR of 12 dB and correlation matrix size of N=60. 

As can be seen from these figures, the algorithm estimates the correct TDOA of 40 
samples or 200 ms for a SNR of 12 dB, and yields a completely incorrect value of -63 
samples for an SNR of 5 dB. 
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I 

Figure 12. Subspace TDOA estimation for an exponentially decaying sinusoid 
using an SNR of 5 dB and correlation matrix size of N=60. 

The rationale for using a large correlation matrix is twofold. First, if the 
correlation matrix is too small, the TDOA estimate tends to be less accurate. A typical 
result is shown in Figure 13 for the 12 dB-SNR case using a correlation matrix of size N 



32 



= 5. As can be seen in this figure, a TDOA of 38 samples (instead of the correct value 40 
samples) was estimated, and the peak is broader than that shown in Figure 11. Secondly, 
for the subspace method to function as a discriminator, it was found that better results 
were achieved at lower SNR, when the correlation matrix and thus the observation space 
was larger. 
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Figure 13. Subspace TDOA estimation for an exponentially decaying sinusoid 
using an SNR of 12 dB and correlation matrix size of N = 5. 

For the bispectrum linear phase detector, no windowing was applied to the data. 
Windowing is usually applied to obtain smooth estimates of the bispectrum [Ref 7:p. 
126]. Due to the short data lengths used in this thesis, however, it was found that 
windowing did not improve the results and was therefore not used. Typical results for the 
bispectrum linear phase detector are shown in Figure 14 and Figure 15 for the 
exponentially decaying sinusoid using SNR values of 12 dB and 5 dB. In both cases the 
TDOA is indicated by the maximum value of the hologram, which in turn is a single 
sample. The performance is similar to that of the subspace method. At 12 dB the result 
is exactly correct while at 5 dB the method gives completely erroneous results. 
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Figure 14. Bispectrum TDOA estimation for an exponentially decaying transient 

using an SNR of 12 dB. 




Figure 15. Bispectrum TDOA estimation for an exponentially decaying transient 

using an SNR of 5 dB. 

Figure 16 shows typical TDOA results using the classical generalized cross- 
correlation methods (GCC) [Ref 2] at a SNR value of 12dB. ROTH [Ref 13], SCOT [Ref 
14] and PHAT [Ref 15] weighting was used in addition to straight cross-correlation (no 
special weighting was used). As can be seen in Figure 16 the correlation, SCOT and 
PHAT algorithms give strong peaks at the estimated TDOA while ROTH gives a very 
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noisy signal with no dominant strong peak. In general, it was found that these techniques 
gave less accurate TDOA estimates than the bispectral and subspace algorithms, even at 
high SNRs. 
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Figure 16. TDOA estimates using (a) cross correlation, (b) ROTH, (c) SCOT and 

(d) PH AT with an SNR of 12 dB. 

TDOA-estimate performance of the subspace and bispectral algorithms was 
evaluated using simulations of one thousand realizations at a fixed SNR. For each 
realization, a different white noise source was added to the original signal. The output of 
these simulations was used to create a statistical measure of performance for TDOA, 
which we call the probability of correct TDOA Pj. Pt represents the fraction of 1000 
trials for which the algorithm estimates the correct time delay between transients received 
at the two sonabuoys. Pt curves were generated for SNR values of 8 dB to 20 dB for 
both the subspace method and the bispectrum linear phase detector. 

Figure 17 and Figure 18 (a) shows the Pt results of the subspace algorithm for the 
CW pulse. As can be seen from the figures, the subspace algorithm gives excellent 
results for SNR >17 dB, where a Pt of 1 .0 is achieved. 
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Figure 17. Pj versus SNR for CW pulse using the subspace linear phase detector. 

The Pt curves for the exponentially decaying transient and the LFM transient are 
shown in Figure 18 (b) and (c). Once again, a Pj value of 1.0 is achieved for SNR > 
17dB. The Pt for the whale transient (Figure 18 (d)) shows poorer results, achieving a 
maximum Pt of 1.0 at a SNR of 20 dB. This poorer performance of the algorithm may 
be be attributed to two factors. First, the whale transient was the longest of all the 
transients used (being 1 s long) thus filling most of the data observation window. 
Secondly, it is a non-linear function giving it a complicated phase structure. It is also 
important to note that the results in Figure 18 (d) are based on a longer data observation 
length than used for the other TDOA results (see Table 3). 

Figure 18 (a) and Figure 19 show the Pt for the CW pulse using the bispectrum 
linear phase detector. From these figures, it can be seen that a PTof 1.0 is never achieved 
but that a Pt of 0.9 is achieved at a SNR of 15 dB. The remaining plots in Figure 18 (b), 
(c) and (d) also show that a Pt of 1.0 is not achieved for the bispectrum method. For the 
exponentially decaying sinusoid (Figure 18 (b)) a PTof 0.9 is achieved at 1 1 dB while for 
the LFM pulse (see Figure 18 (c)) a Pt of 0.9 is achieved at 12 dB. As in the case of the 
subspace methods, the Pt of the whale transient was obtained using a longer data 
observation window, where a maximum Pt of 0.9 was achieved at an SNR of 16 dB 
(Figure 18(d) ). 
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Figure 18. Combined Pt for subspace and bispectrum linear phase detectors: (a) 
CW pulse, (b) exponential decaying sinusoidal transient, (c) LFM pulse and (d) 

whale transient. 
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Figure 19. Pt for the CW pulse using the bispectrum linear phase detector. 
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Figure 18 gives a summary of the TDOA results for all of the experiments. As 
can be seen in Figure 1 8, the bispectrum linear phase detector has higher Pj values at low 
SNR than the subspace linear phase detector for the exponential, LFM and whale 
transients. However at higher values of SNR, the subspace algorithm gives better Pt 
results (i.e., Pt of 1.0) which is never achieved using the bispectrum linear phase 
detector. 

In comparing the two methods (subspace and bispectrum) one may note that 
although the bispectrum method never achieves a Pt of 1 .0, in three out of the four cases 
it reaches the value of Pt = 0.9 earlier than the subspace method. Table 3 compares these 
results by listing the minimum SNR values for which a Pt of 0.9 is achieved in each of 
the test cases. 





CW Pulse 


Exponential 


LFM Pulse 


Whale 


Subspace 


13dB 


15dB 


14dB 


18dB 


Bispectrum 


15dB 


lldB 


12dB 


16dB 



Table 3. Minimum SNR required to achieve Pj = 0.9. 



2. Discrimination Experiments 

The transient discrimination results will be discussed in a fashion similar to that 
of the TDOA results. The discussion will therefore start by giving a brief introduction of 
how the algorithms were implemented by using the CW and LFM pulses as examples. 
After this, the probability curves for discrimination will be discussed for every 
combination of transients. 

As discussed previously, the subspace discriminator makes use of projections into 
multiple subspaces and compares the maximum values of the indicator functions Eq. 

4.11 with the maximum value of the function (Eq. 4.8). It is found that if the 

transients arriving at the sonabuoys are the same the difference in peak values will be 
large as shown in Figure 20(a). If the difference is found to be small (i.e., < 1) then the 
two transients are different, as shown in Figure 20(b). Thus in making these 
comparisons, it is possible to discriminate between transients. 
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Expanded view: 




Figure 20. Typical estimates, (a) when the signals are the same and (b) when 
the signals are different. 

In the case of the bispectrum discriminator, the phase terms (f^R- (ps (see Eq 3.37), 
which are only present in the case of different transient arrivals, add up constructively or 
destructively thus giving features to the hologram of the bispectral linear phase detector. 
For the case where the transients are the same, the phase terms ^ are zero; therefore 
no extra features are added to the hologram. If the difference in magnitude between these 
features is large compared to the average of the hologram, it is possible to discriminate 
between transients. These characteristics can be clearly seen in Figure 21. Figure 21(a) 
shows the single peak that appears for similar transient arrivals while Figme 21(b) shows 
the multiple peaks and extra features that are present due to the phase terms <f)R and ^ for 
different transient arrivals. 
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Figure 21. Bispectrum linear phase detector. Stem plot of |7’(^')| for (a) similar 
transient arrivals and (b) different transient arrivals. 

The experiments for transient discrimination were as follows. Each algorithm 
was evaluated using simulations of one thousand realizations at a fixed SNR. For each 
realization, a different white noise source was added to the original signal. The output of 
these simulations was used to estimate the probability of correct discrimination (PdO- Pdi 
represents the fraction of 1000 trials for which the algorithm discriminates correctly 
between two transients. Pdi curves were generated for SNR values of 10 dB to 20 dB for 
both the subspace and bispectrum discriminators. The probability of incorrect 
discrimination, Pfa, was computed in a similar manner and also used as a performance 
measure. 

Figure 22 shows the results using the subspace discriminator where the 
exponentially decaying sinusoid is present at sonabuoy 1 and various transients are 
present at sonabuoy 2. For these curves, the vertical axis is labeled “Probability” because 
the Pfa and Pdi are shown on the same plot. In Figure 22, the Pdi curve for the exponential 
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Figure 22. Poi and Pfa (using subspace methods) between exponential decaying 
transient and exponentially decaying transient, CW Pulse, noise, LFM pulse and 
whale transient. 



decaying transient arriving at the first and second sonabuoy is shown. It can be seen 
from this curve that the algorithm discriminates correctly with a Poi of 0.9 or larger for 
SNR values greater than 17 dB. Figure 22 also shows the Pfa for the other transient 
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arrivals at the second sonabuoy (i.e., CW pulse, LFM pulse, or the whale transient) and 
the Pfa when only noise (no transient) is present at the second sonabuoy. These curves 
show that the Pfa is veiy low. These are desirable results since they indicate that the 
algorithm can differentiate between the transients with a high probability of correct 
discrimination and a low probability of false alarm. Results for the other transients are 
summarized in Table 4; the corresponding curves are given in appendix A. It is 
important to note that the results for the whale transient are based on a longer data 
observation interval than the others. 



Transient 

Arriving 


B 

U 

O 

Y 

2 


CW Pulse 


Exponential 


LFM Pulse 


Whale 


Noise 


BUOY 1 














CW Pulse 


Poi >0.35 for 
SNR>19dB 


Pfa< 0.08 


Pfa <0.01 


O 

II 


Pfa <0.01 


Exponential 


Pfa < 0.02 


Poi>0.9 for 
SNR>17dB 


Pfa < 0.002 


II 

o 


Pfa < 0.021 


LFM Pulse 


Pfa < 0.03 


Pfa <0.08 


Poi>0.8 for 

SNR>18dB 


Pfa <0.1 


Pfa <0.04 


Whale 




Pfa <0.2 


Pfa <0.18 


Pfa< 0.18 


Poi>0.8 for 
SNR>18dB 


Pfa <0.19 



Table 4. Poi results using subspace methods. 



For the bispectral discriminator, the threshold Or (see Eq. 3.44) was adjusted by 
using a threshold gain to improve the results. The gain was applied to Eq. 3.44 in the 
following way: 



Qt = 



1=0 



N-2 



(5. 13) 



where G is the threshold gain. Threshold gains of 1, 2, 4 and 6 were used. 

Figure 23 shows the discrimination results for the exponentially decaying 
transient using the bispectral discriminator. Again, the vertical axis is labeled 
“Probability” because both the Pfa and Pd; are shown on the same plot. For Figure 23, the 
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transient arriving at sonabuoy 1 is the exponentially decaying sinusoid. The Pqi for this 
transient is shoAvn in Figure 23. The other curves of Figure 23 show the results of the Pfa 
for the CW pulse, noise, LFM pulse and whale transient respectively. The symbols 
shown in Figure 23 are tied to the transient type and are consistant with what is used in 
Appendix B. The subplots in Figure 23 show the results when different threshold gains 
are used. 
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Figure 23. Poi (using high-order methods) between exponentially decaying 
sinusoid transient, arriving at sonabuoy 1 and the CW pulse, exponentially decaying 
Sinusoid, noise, LFM pulse and whale transients arriving at sonabuoy 2. Threshold 
gain values of (a) 1, (b) 2, (c) 4, (d) 6 were used. 
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As can be seen in Figure 23, both Pfa and the Pd, decrease with increasing 
threshold gain. A design compromise must therefore be made to have an acceptably high 
Poi and an acceptably low Pfa. For this thesis, a threshold gain of 4 was chosen. Table 5 
shows the discrimination results for the bispectral discriminator using this value for the 
threshold gain. 



Transient 

Arriving 


B 

U 

O 

Y 

2 


CW Pulse 


Exponential 


LFM Pulse 


Whale 


Noise 


BUOY 1 














CW Pulse 


^Di ^0.8 for 
SNR>17dB 


Pfa <0.2 


Pf3<0.01 


Pfa=0 


Pfa <0.001 


Exponential 


Pfa <0.3 


Poi >0.98 for 
SNR>14dB 


Pfa <0.15 


Pfa = 0 


Pfa <0.001 


LFM Pulse 


Pfa <0.05 


Pfa <0.03 


Poi >0.95 for 
SNR>16dB 


Pfa = 0 


Pfa <0.001 


Whale 




Pfa =0 


Pfa=0 


Pfa = 0 


Poi >0.2 for 
SNR>20dB 


Pfa = 0 



Table 5. Pdi using Bispectrum using a threshold gain of 4. 



The poor discrimination performance of the whale transient given in Table 5 may 
be attributed to the relatively short observation time of the data. Better results may have 
been obtained if longer observation times were used; however, it was not practical to run 
a large number of simulations for these longer observation times. 

C. SUMMARY OF RESULTS 

The TDOA and discriminator results are summarized in Table 6 including the 
advantages and disadvantages of both algorithms. As has been previously shown in 
Figure 1 8, the bispectral linear phase detector tends to give a larger number of correct 
TDOA estimates at low SNR. On the other hand, the subspace algorithm gives 
consistently correct results (Pj = 1.0) at high SNR, which is never achieved by the 
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bispectrum algorithm for SNR in the range of 8 dB to 20 dB. Therefore the more 
desirable result may depend on the application where the algorithms will be used. 

For discrimination, the bispectral method gave the highest Pdi results when a 
threshold gain of 4 was used. However, the Pfa increased with increasing SNR (see 
Figure 23). This is an undesirable result since even transients with high SNR will not be 
able to be separated. On the other hand, the subspace discriminator gave lower Poi values 
but had the advantage that the Pfa remained constant for all SNR values tested. 



Algorithm 


TDOA 


Discrimination 


Bispectrum 


Best Result: SNR > lldB 
gives a Pt of 0.9 
(exponential transient) 
Disadvantage: A Pt of 1 is 
never achieved at SNRs < 
20dB 

Advantage: High Pj’s are 
achieved at low SNRs 


Best Result: SNR > 14dB gives a 
Poi of 0.98 (exponential transient). 
Disadvantages: Results depend on 
choice of threshold gain. Pfa 

increases with increasing SNR. 
Advantages: Can achieve a high 
Poi at relatively low SNRs. Has a 
low Pfa for noise 


Subspace 


Best Result: SNR > 17dB 
gives a Pt of 1. (exponential 
transient) 

Disadvantage: High Pt’s 

are only achieved at high 
SNRs 

Advantage: Pt’s of 1 are 

achieved. 


Best Result: SNR > 17dB gives a 
Poi of 0.9 (exponential transient). 
Disadvantages: Computationally 

intensive. Low Po, at high SNRs 
Advantages: Low Pfa for all 

transients 



Table 6. Summary of TDOA and discrimination results. 
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VI. CONCLUSIONS 



A. THESIS SUMMARY 

In this thesis, we have developed and compared two algorithms, namely the 
bispectrum and subspace linear phase detectors. These algorithms were developed for 
the purposes of transient discrimination and TDOA estimation, in order to be used as part 
of a transient tool suite to aid in the estimation of a submarine’s position. Chapters III 
and IV of this thesis provide an analysis of these two algorithms. Two performance 
measures were used to evaluate the algorithms, namely the probability of correct TDOA 
(Pt) and the probability of correct discrimination (Poi)- 

Chapter V discusses the simulations and transients used in the evaluation of the 
algorithms. Of the four transients used, it was found that the whale transient, which is the 
longest transient, gave the worst performance. This may be due to the relatively short 
data observation times used in evaluating the algorithms. It was found that longer 
observation times produced better results for this transient; however, it was impractical to 
run a large number of cases at the longer observation times. 

It was found that both algorithms had advantages and disadvantages as 
summarized in Table 6. In general, it can be said that for TDOA the bispectral linear 
phase detector gave better results at low SNR while the subspace linear phase detector 
worked better at the higher SNR. 

For discrimination, it was found that the bispectral discriminator gave higher Poi 
than the subspace discriminator. However, the Pfa of the bispectral discriminator 
increased at higher SNR while the subspace discriminator gave a constant Pfa for all SNR 
values tested. For discrimination, the advantage of having a constant Pfa is desirable. 
Therefore the subspace discriminator is the best option even although it produced lower 
Poi than the bispectral discriminator. As discussed in Chapter V, there are design trade 
offs between processing speed and performance that need to be made. For the bispectral 
linear phase detector, this trade off is in terms of threshold gain; for the subspace linear 
phase detector, this trade off is in terms of correlation matrix size. 
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B. 



FUTURE WORK 



The results of this thesis lead to some interesting topics for future research. These 
topics include; 

(1) The effects of the environment on a transient . In the evaluation of the two 
algorithms conducted in this thesis, no work was done to determine the sensitivity of the 
algorithms to changes in transient shape due to the environment. To do this, a 
propagation model must be used. 

(2) The effects of colored noise . In the simulations additive white Gaussian noise 
was used. This is a poor reflection of reality since ocean noise is frequently colored [Ref 
22]. It will therefore be important to see how these algorithms perform when additive 
colored noise is used. 

(3) Operational testing . The algorithms were tested using synthetic data. 

Ultimately, there is a need to test them on real data to obtain a measure of their 
performance and applicability in the field. 
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APPENDIX A 



Pdi and Pfa using SUBSPACE METHOD 
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Figure 24. Pdi for the CW transient arriving at the first sonabuoy. 
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Figure 25. 



Poi for the LFM transient arriving at the first sonabuoy. 
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Figure 26 . Pdi for the whale transient arriving at the first sonabuoy. 
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APPENDIX B 



Pdi and Pfa using BISPECTRUM METHOD 
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Figure 27. Pdi (using bispectral methods) between CW pulse transient , arriving 
at sonabuoy 1 and the CW pulse, exponentially decaying sinusoid, noise, LFM pulse 
and whale transients arriving at sonabuoy 2. Threshold gain values of (a) 1, (b) 2, 
(c) 4, (d) 6 were used. 
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Figure 28. Pd, (using bispectral methods) between LFM pulse transient , arriving 
at sonabuoy 1 and the CW pulse, exponentially decaying sinusoid, noise, LFM pulse 
and whale transients arriving at sonabuoy 2. Threshold gain values of (a) 1, (b) 2, 
(c) 4, (d) 6 were used. 
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Figure 29. Ppi (using bispectral methods) between whale transient , arriving at 
sonabuoy 1 and the CW pulse, exponentially decaying sinusoid, noise, LFM pulse 
and whale transients arriving at sonabuoy 2. Threshold gain values of (a) 1, (b) 2, 
(c) 4, (d) 6 were used. 
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